package ij3d.utils;

import ij.*;
import ij.measure.CurveFitter;
import ij.process.*;

/**
 *  class allowing : <BR>
 *  storage of arrays with util function such as statistics
 *  and sorting
 *
 *@author     Cedric MESSAOUDI & Thomas BOUDIER
 *@created    17 septembre 2003
 */
public class ArrayUtil {

    private float tab[];
    private int size;
    private boolean sorted;
    private boolean debug = false;

    /**
     *  constructeur
     *
     *@param  size  number of elements
     */
    public ArrayUtil(int size) {
        this.size = size;
        tab = new float[size];
        for (int i = 0; i < size; i++) {
            tab[i] = 0f;
        }
        sorted = false;
    }

    /**
     *  constructeur
     *
     *@param  data  float array
     */
    public ArrayUtil(float[] data) {
        this.size = data.length;
        sorted = false;
        tab = data;
    }

    /**
     *  put a value to a index
     *
     *@param  pos  position in the array
     *@param  value     value to put
     *@return          false if position does not exist
     */
    public boolean putValue(int pos, float value) {
        if (pos < size) {
            tab[pos] = value;
            sorted = false;
            return true;
        } else {
            return false;
        }
    }

    /**
     *  get value at an index
     *
     *@param  pos  position in the array
     *@return          the value, or NaN if does not exist
     */
    public float getValue(int pos) {
        if (pos < size) {
            return tab[pos];
        } else {
            return Float.NaN;
        }
    }

    /**
     *  get the number of elements
     *
     *@return   number of elements
     */
    public int getSize() {
        return size;
    }

    /**
     *  Gets the array attribute of the ArrayUtil object
     *
     *@return    The array value
     */
    public float[] getArray() {
        return tab;
    }

    /**
     * new size of the array (can incresase size of array)
     *
     *@param  size  The new size value
     */
    public void setSize(int size) {
        if (size > this.size) {
            float[] temp = new float[size];
            for (int i = 0; i < this.size; i++) {
                temp[i] = tab[i];
            }

            tab = temp;
        }
        this.size = size;
    }

    /**
     *  The maximum value
     *
     *@return   max value
     */
    public float getMaximum() {
        float max = tab[0];
        for (int i = 1; i < size; i++) {
            if (tab[i] > max) {
                max = tab[i];
            }
        }

        return max;
    }

    /**
     * Last index whose value is not zero
     *
     *@return    The limitSup value
     */
    public int getLimitSup() {
        int i = size - 1;
        while ((i >= 0) && (tab[i] == 0)) {
            i--;
        }

        return i;
    }

    /**
     * First encoutered value not zero, starting from th
     *
     *@param  th  minimum index
     *@return     The limitInf value
     */
    public int getLimitInf(int th) {
        int i = th;
        while ((i < size) && (tab[i] == 0)) {
            i++;
        }

        return i;
    }

    /**
     *  The minimum value
     *
     *@return   min value
     */
    public float getMinimum() {
        float min = tab[0];
        for (int i = 1; i < size; i++) {
            if (tab[i] < min) {
                min = tab[i];
            }
        }

        return min;
    }

    /**
     *  Gets the minimum value above a threshold
     *
     *@param  th  the threshold
     *@return     The minimumAbove value
     */
    public float getMinimumAbove(float th) {
        float min = Float.MAX_VALUE;
        for (int i = 0; i < size; i++) {
            if (tab[i] > th) {
                if (tab[i] < min) {
                    min = tab[i];
                }
            }
        }
        if (min < Float.MAX_VALUE) {
            return min;
        } else {
            return th;
        }
    }

    /**
     *  Average value
     *
     *@return   average value
     */
    public float getMean() {
        float total = 0;
        for (int i = 0; i < size; i++) {
            total += tab[i];
        }

        return total / size;
    }

    /**
     * sum of all elements
     *
     *@return   sum
     */
    public float getSum() {
        float total = 0;
        for (int i = 0; i < size; i++) {
            total += tab[i];
        }

        return total;
    }

    /**
     *  standard deviation
     *
     *@return   std dev
     */
    public float getSigma() {
        float var = getVariance();
        return (float) Math.sqrt(var);
    }

    /**
     * skewness, order 3
     *
     * @return skewness value (- for right distribution, 0 symmetric, + left distribution)
     */
    public float getSkewness() {
        float mu = getMean();
        float si = getSigma();
        float diff;

        float total = 0;
        for (int i = 0; i < size; i++) {
            diff = (tab[i] - mu) / si;
            total += diff * diff * diff;
            //System.out.println(tab[i]);
        }

        if (size > 1) {
            {
                float fsize = (float) size;
                total *= fsize / ((fsize - 1) * (fsize - 2));
            }
        } else {
            total = (float) 0.0;
        }
        return (float) total;
    }

    /**
     * Kurtosis, order 4
     *
     * @return kurtosis (peakiness)
     */
    public float getKurtosis() {
        float mu = getMean();
        float var = getVariance();
        float diff;

        float total = 0;
        for (int i = 0; i < size; i++) {
            diff = tab[i] - mu;
            total += diff * diff * diff * diff;
        }

        if (size > 1) {
            float fsize = (float) size;
            total *= ((fsize) * (fsize + 1)) / ((fsize - 1) * (fsize - 2) * (fsize - 3) * var * var);
            total -= 3 * (fsize - 1) * (fsize - 1) / ((fsize - 2) * (fsize - 3));

        } else {
            total = 0;
        }

        return (float) total;
    }

    /**
     * Variance value (not biased)
     *
     * @return variance
     */
    public float getVariance() {
        if(size==1) return 0;

        float mean = this.getMean();
        float total = 0;
        for (int i = 0; i < size; i++) {
            total += (tab[i] - mean) * (tab[i] - mean);
        }

        if (size > 1) {
            total /= (size - 1);
        } else {
            total = 0;
        }

        return (float) total;
    }

    /**
     * Variance value (not biased)
     * fast version
     *
     * @return variance
     */
    public float getVariance2() {
        if(size==1) return 0;
        
        double total = 0;
        double total2 = 0;

        for (int i = 0; i < size; i++) {
            total += tab[i];
            total2 += tab[i] * tab[i];
        }

        float var = (float) ((total2 - (total * total / size)) / (size - 1));
        return var;
    }

    /**
     * Array with only distinct values (sorted)
     *
     * @return
     */
    public ArrayUtil distinctValues() {
        this.sort();
        ArrayUtil V = new ArrayUtil(this.getSize());
        int s = 0;
        float tmp = this.getValue(0);
        V.addValue(0, tmp);
        s++;
        int p = 1;
        int si = this.getSize();
        while (p < si) {
            while ((p < si) && (this.getValue(p) == tmp)) {
                p++;
            }
            if (p < si) {
                tmp = this.getValue(p);
                V.addValue(s, tmp);
                s++;
                p++;
            }
        }

        return V.getSubTabUtil(0, s);
    }

    /**
     * local mean filtering
     *
     *@param  rad  radius for filtering
     *@return      filtered array
     */
    public ArrayUtil localMean(int rad) {
        ArrayUtil tab2 = new ArrayUtil(size);
        for (int i = 0; i < size; i++) {
            if (i > rad && i < size - rad) {
                int total = 0;
                int nombre = 0;
                for (int j = i - rad; j <= i + rad; j++) {
                    total += tab[j];
                    nombre++;
                }
                tab2.putValue(i, total / nombre);
            } else {
                tab2.putValue(i, tab[i]);
            }
        }

        return tab2;
    }

    /**
     *  Adds a value at a index
     *
     *@param  position  the position index
     *@param  value     the value to add
     *@return           true if position < size
     */
    public boolean addValue(int position, float value) {
        if (position < size) {
            tab[position] += value;
            sorted = false;
            return true;
        }
        return false;
    }

    /**
     * Sort the array (Shell-Meitzner algorithm)
     */
    public void sort() {
        float aux;
        for (int ecart = size / 2; ecart > 0; ecart /= 2) {
            int i = 0;
            int pos_mem = 0;
            while (i < (size - ecart)) {
                if (tab[i] > tab[i + ecart]) {
                    aux = tab[i];
                    tab[i] = tab[i + ecart];
                    tab[i + ecart] = aux;
                    if (i - ecart >= 0) {
                        i -= ecart;
                    } else {
                        i = ++pos_mem;
                    }

                } else {
                    i = ++pos_mem;
                }
            }

        }
    }

    /**
     *  Get the indexes of sorted array (Shell-Meitzner algorithm)
     * do not sort the array
     *
     *@return    the sorted index
     */
    public int[] sortIndex() {
        int td[] = new int[size];
        for (int i = 0; i < size; i++) {
            td[i] = i;
        }
        int aux;
        for (int ecart = size / 2; ecart > 0; ecart /= 2) {
            int i = 0;
            int pos_mem = 0;
            while (i < (size - ecart)) {
                if (tab[td[i]] > tab[td[i + ecart]]) {
                    aux = td[i];
                    td[i] = td[i + ecart];
                    td[i + ecart] = aux;
                    if (i - ecart >= 0) {
                        i -= ecart;
                    } else {
                        i = ++pos_mem;
                    }

                } else {
                    i = ++pos_mem;
                }
            }
        }
        return td;
    }

    /**
     *  The median (sorted array)
     *
     *@return   mediane
     */
    public float medianSort() {
        if (!sorted) {
            sort();
        }

        return tab[size / 2];
    }

    /**
     *  Find median. Modified algorithm according to
     *  http://www.geocities.com/zabrodskyvlada/3alg.html Contributed by HeinzKlar.
     *  (copied form ij.plugin.filter.rankfilters
     *
     *@return    Median value
     */
    public float median() {
        final int nValues = size;
        final int nv1b2 = (nValues - 1) / 2;
        int i;
        int j;
        int l = 0;
        int m = nValues - 1;
        float med = tab[nv1b2];
        float dum;

        while (l < m) {
            i = l;
            j = m;
            do {
                while (tab[i] < med) {
                    i++;
                }
                while (med < tab[j]) {
                    j--;
                }
                dum = tab[j];
                tab[j] = tab[i];
                tab[i] = dum;
                i++;
                j--;
            } while ((j >= nv1b2) && (i <= nv1b2));
            if (j < nv1b2) {
                l = i;
            }
            if (nv1b2 < i) {
                m = j;
            }
            med = tab[nv1b2];
        }
        return med;
    }

    /**
     *  Convolution with a kernel
     *
     *@param  kernel  the kernel
     *@param  D       the divisor
     *@return         the convolved value
     */
    public float convolve(float[] kernel, float D) {
        float sum;
        if (size != kernel.length) {
            return -1;
        } else {
            sum = 0;
            for (int i = 0; i < size; i++) {
                sum += tab[i] * kernel[i];
            }
            sum /= D;
        }
        return sum;
    }

    /**
     * Index to get x % of the values
     *
     *@param  percent  percentage(1.0=100; 0.5=50% ...)
     *@return        index in the array
     */
    public int indexOfSumPercent(double percent) {
        int total = 0;
        int sum = 0;
        for (int i = 0; i < size; i++) {
            total += tab[i];
        }

        for (int i = 0; i < size; i++) {
            sum += tab[i];
            if ((double) sum >= percent * total) {
                return (i > 0 ? i - 1 : 0);
            }
        }
        return size - 1;
    }

    /**
     *  The value is in the array ?
     *
     *@param  val  the value
     *@return      true if value in the array
     */
    public boolean hasValue(int val) {
        int i = 0;
        boolean in = false;
        while ((!in) && (i < size)) {
            if (tab[i] == val) {
                in = true;
            }

            i++;
        }
        return in;
    }

    /**
     *  test if array has only a particular value
     *
     *@param  val  the value
     *@return      true if only this value is inside the array
     */
    public boolean hasOnlyValue(int val) {
        int i = 0;
        while (i < size) {
            if (tab[i] != val) {
                return false;
            }
            i++;
        }
        return true;
    }

    /**
     *  isodata algorithm (for histogram)
     *
     *@return    the computed isodata value
     */
    public int IsoData() {
        int sup = this.getLimitSup();
        int m1 = sup / 2;
        int m0 = 0;
        float md;
        float mg;
        float nombre_pixels;
        int nb = 0;

        while ((Math.abs(m0 - m1) > 2) && (nb < size)) {
            nb++;
            m0 = m1;
            md = 0;
            nombre_pixels = 0;
            if (debug) {
                IJ.log("size=" + size + " m0=" + m0);
            }

            for (int i = 0; i < m0; i++) {
                nombre_pixels = nombre_pixels + tab[i];
                md = md + tab[i] * i;
            }
            if (debug) {
                IJ.log("md=" + md + " nb=" + nombre_pixels);
            }

            md = md / nombre_pixels;
            mg = 0;
            nombre_pixels = 0;
            for (int i = m0; i <= sup; i++) {
                nombre_pixels = nombre_pixels + tab[i];
                mg = mg + tab[i] * i;
            }
            if (debug) {
                IJ.log("mg=" + mg + " nb=" + nombre_pixels);
            }

            mg = mg / nombre_pixels;
            m1 = (int) ((md + mg) / 2);
            if (debug) {
                IJ.log("m0=" + m0 + " m1=" + m1);
            }
        }
        if (nb == size) {
            if (debug) {
                System.out.println("Pb convergence moyenne intermediaire");
            }
        }
        return m1;
    }

    /**
     *  drawing of the array as an image
     *
     *@param  scale  scale value
     */
    public void draw(int scale) {
        int haut = 256;
        float val = 0;
        float max = getMaximum();
        ByteProcessor ip = new ByteProcessor(size, haut);
        ip.setColor(255);
        for (int i = 0; i < size; i++) {
            ip.moveTo(i, haut - 1);
            if (scale == -1) {
                val = (haut * tab[i] / max);
            } else if (tab[i] > scale) {
                val = haut - 1;
            } else {
                val = (haut * tab[i] / scale);
            }

            ip.lineTo(i, (int) (haut - val));
        }
        ImagePlus iplus = new ImagePlus("tab", ip);
        iplus.show();
    }

    /**
     *  random mix of the values
     */
    public void randomize() {
        sorted = false;
        for (int i = 0; i < size - 1; i++) {
            int pos = (int) Math.round(Math.random() * (size - 1 - i) + i);
            float aux = tab[i];
            tab[i] = tab[pos];
            tab[pos] = aux;
        }
    }

    /**
     *  put values inside an array
     *
     *@param  tmp  the array of values to put
     */
    public void putValues(ArrayUtil tmp) {
        if (tmp.getSize() > tab.length) {
            tab = new float[tmp.getSize()];
        }

        this.setSize(tmp.getSize());
        for (int i = 0; i < size; i++) {
            tab[i] = tmp.getValue(i);
        }
    }

    /**
     *  test if two arrays are identical
     *
     *@param  other  the other array
     *@return        true if arrays are equal
     */
    public boolean isEqual(ArrayUtil other) {
        if (other.getSize() != this.getSize()) {
            return false;
        }
        return isEqual(other, 0, this.getSize() - 1);
    }

    /**
     *  test if two arrays are identical between two indexes
     *
     *@param  other  the other array
     *@param  begin  start index
     *@param  end    end index
     *@return        The equal value
     */
    public boolean isEqual(ArrayUtil other, int begin, int end) {
        boolean result = true;
        for (int i = begin; i <= end; i++) {
            if (this.tab[i] != other.tab[i]) {
                result = false;
            }
        }

        return result;
    }

    /**
     *  Gets the copy of the ArrayUtil object
     *
     *@return    The copy array
     */
    public ArrayUtil getCopy() {
        ArrayUtil result = new ArrayUtil(this.size);
        for (int i = 0; i < this.size; i++) {
            result.tab[i] = this.tab[i];
        }
        return result;
    }

    /**
     *  information to be displayed
     *
     *@return    text
     */
    @Override
    public String toString() {
        String str = "{" + tab[0];
        for (int i = 1; i < size; i++) {
            str = str + ", " + tab[i];
        }
        return str + "}";
    }

    /**
     *  remove a value at a position
     *
     *@param  index  the position index
     */
    public void removeValueAt(int index) {
        size--;
        for (int i = index; i < size; i++) {
            tab[i] = tab[i + 1];
        }
    }

    /**
     *  concat an array to the end
     *
     *@param  tabToAdd  the array to concat
     */
    public void concat(ArrayUtil tabToAdd) {
        int newsize = size + tabToAdd.size;
        float[] tmp = new float[newsize];
        for (int i = 0; i < size; i++) {
            tmp[i] = tab[i];
        }
        for (int i = 0; i < tabToAdd.size; i++) {
            tmp[size + i] = tabToAdd.tab[i];
        }
        this.tab = tmp;
        size = newsize;
    }

    /**
     *  extract a sub-array from this array
     *
     *@param  startindex  start index
     *@param  newsize     new size
     *@return             The subTabUtil value
     */
    public ArrayUtil getSubTabUtil(int startindex, int newsize) {
        ArrayUtil tmp = new ArrayUtil(newsize);
        for (int i = 0; i < newsize; i++, startindex++) {
            tmp.putValue(i, getValue(startindex));
        }
        return tmp;
    }

    /**
     * Fit a gauusian to the values (radial distribution)
     *
     * @param values to fit
     * @param maxR maximum interval
     * @return the gaussian fitted values
     */
    public static double[] fitGaussian(double[] values, int maxR) {
        // create values
        double[] id = new double[values.length];
        int c = 0;
        double minVal = values[0];
        double maxVal = values[0];
        double val;
        for (int i = -maxR; i <= maxR; i++) {
            val = values[c];
            if (val > maxVal) {
                maxVal = val;
            }
            if (val < minVal) {
                minVal = val;
            }
            id[c] = i;
            c++;
        }
        IJ.log("min=" + minVal + " max=" + maxVal);
        // fitting by a gaussian
        CurveFitter fit = new CurveFitter(id, values);
        double[] params = new double[4];
        params[0] = minVal;
        params[1] = maxVal;
        params[2] = 0;
        params[3] = 3;
        IJ.log("res=" + params[0] + " " + params[1] + " " + params[2] + " " + params[3]);
        fit.setInitialParameters(params);
        fit.setMaxIterations(10000);
        fit.setRestarts(1000);
        fit.doFit(CurveFitter.GAUSSIAN);
        params = fit.getParams();
        IJ.log("res=" + params[0] + " " + params[1] + " " + params[2] + " " + params[3]);
        IJ.log("error=" + fit.getRSquared());

        return params;
    }
}
